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New Measurements of the Radio Photosphere of Mira based on 
Data from the JVLA and ALMA 

L. D. Matthews^ M. J. Reid^, & K. M. Menten^ 

ABSTRACT 

We present new measurements of the millimeter wavelength continuum emis¬ 
sion from the long period variable Mira (o Ceti) at frequencies of 46 GHz, 96 GHz, 
and 229 GHz (A ~7 mm, 3 mm, and 1 mm) based on observations obtained 
with the Jansky Very Large Array (JVLA) and the Atacama Large Millime¬ 
ter/submillimeter Array (ALMA). The measured millimeter flux densities are 
consistent with a radio photosphere model derived from previous observations, 
where flux density Sy oc The stellar disk is resolved, and the measurements 
indicate a decrease in the size of the radio photosphere at higher frequencies, as 
expected if the opacity decreases at shorter wavelengths. The shape of the radio 
photosphere is found to be slightly elongated, with a flattening of ~10-20%. The 
data also reveal evidence for brightness non-uniformities on the surface of Mira 
at radio wavelengths. Mira’s hot companion, Mira B was detected at all three 
observed wavelengths, and we measure a radius for its radio-emitting surface of 
2.0 X 10^® cm. The data presented here highlight the power of the JVLA and 
ALMA for the study of the atmospheres of evolved stars. 

Subject headings: stars: AGB and post-AGB - stars: atmospheres - stars: fun¬ 
damental parameters - stars: imaging 


1. Introduction 

Mira (o Geti) is the archetype of its class of long-period variable stars and is one of 
the most extensively studied stars in the Galaxy. Mira has an oxygen-rich chemistry and 
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is undergoing mass loss at a rate of ~ 2 X 10 ^ Mq yr ^ (Ryde & Schoier 2001)ll| The 
pulsation period of Mira is 332 days, during which the visible light of the star varies by up 
to ~8 magnitudes. Mira is also a member of the nearest symbiotic binary system, Mira AB. 
The companion, Mira B, has a projected separation of ~0"5 (55 AU; Karovska et ah 1997; 
Ramstedt et ah 2014) and is likely to be a white dwarf (Sokoloski & Bildsten 2010). 

Both components of the Mira AB binary are detectable at radio wavelengths (Matthews 
& Karovska 2006). Matthews & Karovska showed that the radio emission from Mira B is 
consistent with free-free emission originating from an ionized circumstellar region ~ 10^^ — 
10^^ cm in radius. Based on measurements of the flux density Sb,u between 8.5 and 43.3 GHz, 
they derived a £t to the radio spectrum of the form Sb,i> = (0.010 ± 0.008)z/q^®'^°'^® mJy. 

The radio emission from Mira A has an origin quite different from its companion. Reid 
& Menten (1997; hereafter RM97) showed that the centimeter wavelength emission detected 
from Mira A and other nearby, long-period variables can be explained by the existence of 
a “radio photosphere” with a radius approximately twice that of the classical photospheric 
radius i?* (dehned by the line-free regions of the optical-infrared spectrum). For Mira-type 
variables, R* is typically ~l-2 AU, while between ~1-2R* is a molecular layer (Perrin et ah 
2004) that may be nearly opaque in the visible and infrared. The radio photosphere thus 
lies adjacent to this molecular layer. 

The emission arising from the radio photosphere has a spectral index a ~2 (consis¬ 
tent with optically thick blackbody emission) all the way from centimeter to submillimeter 
wavelengths. However, the flux densities over this spectral range are roughly a factor of two 
higher than expected from an extrapolation of the optical photospheric emission. Based on 
the model developed by RM97, the dominant source of opacity in the radio photosphere 
is the interaction between free electrons (mostly from ionization of K and Na) and neutral 
H and H 2 . At millimeter and centimeter wavelengths, unity optical depth is achieved at a 
radius of ~3 to 4 AU, where the density and temperature are 10^^ cm ^ and ~ 1600 K, 
respectively. 

The detection and characterization of radio photospheres in asymptotic giant branch 
(AGB) stars provide a new tool for the investigation of the critical region between the 
classic photosphere and the dust formation zone where the stellar wind is launched. In 
oxygen-rich AGB stars, the radio photosphere also lies just interior to the atmospheric region 
where SiO masers arise, and contemporaneous observations of these maser lines provide 
additional diagnostics of the density, temperature, dynamics, and shock propagation within 


^Throughout the paper we adopt a distance of 110±9 pc (Haniff et al. 1995), based on the period- 
luminosity relation of Feast et al. 1989. 
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this complex and important atmospheric region (RM97; Gray et ah 2009). As noted by 
RM97, red giants with well-characterized radio photospheres may also be useful as calibration 
sources in the millimeter and submillimeter regimes. 

The characteristic sizes and temperatures of radio photospheres imply that resolved 
imaging observations are within the reach of current radio interferometers for nearby stars. 
Using the Very Large Array (VLA) of the National Radio Astronomy Observatorjo in its 
highest resolution configuration, Reid & Menten (2007) obtained resolved images of three 
oxygen-rich long-period variables (including Mira) at a wavelength of 7 mm, while Menten et 
ah (2012) resolved the carbon star IRC-l-10216. One intriguing finding from these studies is 
evidence that the radio photospheres in some cases appear distinctly non-circular. However, 
since observations were available at only a single epoch and at a single point during the stellar 
pulsation cycle, it was difficult to distinguish between possible origins for these observed 
shapes, and it is unknown whether the shapes are static or time-variable. 

The exquisite continuum sensitivity of the recently upgraded Jansky VLA (JVLA; Perley 
et al. 2011) and the Atacama Large Millimeter/submillimeter Array (ALMA; Wootten & 
Thompson 2009) now provide a tremendous leap in our ability to obtain high-quality images 
of radio photospheres over a broad range of wavelengths. At 7 mm, the continuum sensitivity 
of the JVLA exceeds that of the pre-upgrade VLA by an order of magnitude. At the same 
time ALMA supplies sub-mJy sensitivity with only modest integration times, and its longest 
baselines are capable of providing angular resolution comparable to, or higher than the JVLA, 
depending on observing frequency. Here we present an analysis of new observations of Mira 
obtained with the JVLA and ALMA at wavelengths of 7 mm (46 GHz) 3 mm (96 GHz), and 
1 mm (229 GHz). The quality of the data underscore the potential of these two instruments 
for the study of stellar atmospheres. All of the data were obtained during the course of a 
nine-month period (i.e., during a timescale less than a single stellar pulsation cycle) in 2014. 


2. Observations 
2.1. JVLA Data 

Observations of Mira at 46 GHz (7 mm) were carried out with the JVLA in its most 
extended (A) configuration (0.68 km to 36.4 km baselines) on 2014 February 23 during 
a 3-hour observing session. Data were recorded with 2-second integration times. Offset 
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corrections as refinements to the reference antenna pointing model were evaluated hourly 
using observations of a strong point source at 8 GHz. 3C48 was observed as an absolute 
flux calibrator, and the bright continuum source J0359+5057 (NRAO150) was observed as 
a bandpass calibrator. 

To obtain maximum continuum sensitivity, a 3-bit observing mode was employed, and 
the WIDAR correlator was conhgured with 4 baseband pairs tuned to contiguously cover the 
frequency range from 41.8-49.4 GHz with dual circular polarizations. Each baseband pair 
contained 16 subbands, each with a bandwidth of 128 MHz and 128 spectral channels. This 
frequency range included the SiO v=l and n=2, J=l-0 maser transitions (see below). The 
total on-source integration time for Mira was approximately 1 hour 51 minutes, split into 
two observing blocks that were bracketed by observations of the gain (amplitude and phase) 
calibrator J0215-0222. 


2.2. ALMA Data 

ALMA observations of Mira were performed in Band 3 (94 GHz) on 2014 October 17 and 
25 and in Band 6 (229 GHz) on 2014 October 29 and 2014 November 1, respectively using 
an array conhguration with baselines up to 15.2 km. The data were obtained by ALMA staff 
for science verihcation purposes (ALMA Partnership 2015) and were released to the public 
in 2015 February. A total of 38 (39) antennas were present in the array for the Band 3 (6) 
observations. The total on-source integration times for Mira were approximately 1.6 hours 
in Band 3 and 1.1 hours in Band 6. Data from the different observing dates were combined 
for each band in the present analysis. 


3. Data Processing 
3.1. JVLA Data 


Data processing for the JVLA data was performed using the Astronomical Image Pro¬ 
cessing System (AIPS; Greisen 2003), taking into account special considerations necessary 
for handling JVLA data as outlined in Appendix E of the AIPS Gookbookj^. 

The data were loaded into AIPS directly from archival science data model (ASDM) 
format hies using the Obit software package (Gotton 2008). The default calibration (‘GL’) 


^http: //www.aips.nrao.edu/cook.html 
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table was then recreated to update the single-dish antenna gain curve and atmospheric 
opacity corrections. 

After updating the antenna positions and flagging obviously corrupted data, a requan¬ 
tizer gain correction was applied, as required for 3-bit data. Next, a fringe £t was per¬ 
formed using one minute of data from the observation of J0359-I-5057 to correct the in¬ 
strumental delay settings. Separate delay solutions were determined for each of the four 
independent basebands. Bandpass calibration was performed using NRAO150, assuming 
a spectral index a=0.4 (Trippe et al. 2010). The absolute flux density scale was cal¬ 
culated using the latest time-dependent flux density values for 3C48 from Perley & Butler 
(2013). The resulting flux density Siy as a function of frequency for 3C48 is therefore given by 
log(S'i.) = 1.3324 —0.7690(log(z/)) — 0.1950(log(i^))^-1-0.059(log(z/))^, where u is the frequency 
expressed in GHz. The frequency-independent portion of the complex gains (amplitude and 
phase) was calibrated using J0215-0222 to remove possible slow (hour timescale) instrumen¬ 
tal gain drifts. Additionally, a correction of —3.324" in DEC was applied to the observed 
(J2000.0) position of Mira to account for the large proper motion of the star. This was re¬ 
quired after the discovery that corrections for parallax and proper motion were inadvertently 
excluded from the source position entered in the observation scheduling hie. 

Mira is a well-known source of SiO maser emission (e.g.. Cotton et al. 2006; McIntosh 
& Bonde 2010). We therefore employed the method developed by Reid & Menten (1990) to 
improve the phase calibration using the stellar maser signal. After completion of the initial 
calibration steps described above, the spectral channel containing the brightest SiO n=l, 
J=l-0 maser emission was split from the main data set, and several iterations of phase-only 
self-calibration were performed on the maser data until convergence was reached. From 
these solutions, phase corrections appropriate for each of the 64 subbands comprising the 
full 8 GHz bandwidth of the continuum measurements were derived using the AIPS task 
SNP2D. (Effectively, the phase correction 0j,/ for a given subband with central frequency u' 
then becomes = {v'/vo)(t)o, where z/q is the frequency of the reference channel and 0o is 
the computed phase correction in that channel). If the maser emission is sufficiently strong, 
in principle this approach allows achievement of nearly perfect “seeing” in the resulting 
continuum measurements (Reid & Menten 1990; RM97). However, at the phase Mira was 
observed, the SiO masers were in a relatively low state (integrated hux density 7 Jy), 
and the modest signal-to-noise ratio on some of the longest baselines was a limiting factor 
in the calibration solutions (see Appendix). 

After applying the derived phase corrections, a second round of self-calibration was 
performed on the reference channel, this time solving for both the amplitudes and the phases. 
To mitigate possible drift in the amplitudes, the gains were normalized during this step, and 
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the resulting change in the amplitude scale was measured to be less than 0.5%. Once these 
gain corrections had been applied to the full data set, spectral channels 1-3 and 125-128 of 
each subband (where noise is signihcantly higher) were flagged, along with spectral channels 
containing line emission. Next, the data were averaged to a time resolution of 10 seconds 
to reduce the data volume and Hanning smoothed to reduce ringing caused by the narrow 
maser lines. At this stage, optimized weights for the visibility data were computed using the 
AIPS task REWAY, and hnally, the data were further averaged in frequency, resulting in 8 
spectral channels in each of the 64 subbands (IFs). 


3.2. ALMA Data 

We retrieved the pipeline calibrated u-v data in the form of CASA measurement sets 
for each of the two ALMA bands from the ALMA Science Verihcation Data web pages0 
Details of the processing steps used to produce these data sets are available at that site. 

The flux density scale for the Band 3 data taken on 17 October was calibrated using 
the source J0334-4008, assuming a flux density at 86.2251 GHz of 1.6589 Jy and a spectral 
index a = —0.7146. For the Band 3 data from 25 October, J0238-I-1636 and J2258-2758 
were used as flux density calibrators, assuming respectively Si, =1.8550 Jy, a = —0.269071 
and S'i,=1.1315 Jy, a = —0.7282. For Band 6, the flux density scale for the 29 October 
data was set using J0334-4008, assuming S'y=0.8282 Jy at 229.5497 GHz and a = —0.7029. 
For the 01 November data, the flux calibrator was J0238-I-1636, adopting S'y=1.3893 Jy and 
a = -0.2005 at 229.5485 GHz. 

At this stage, the continuum portion of the Band 3 data included a total of nine spectral 
windows (three in each of three independent basebands) and dual linear polarizations. Each 
spectral window contained eight spectral channels (250 MHz each) across a ~1.8 GHz band 
(i.e., with slight spectral overlap between channels). Adjacent channels thus have slight 
spectral overlap. The frequency ranges sampled were 87.3-89.0 GHz, 97.4-99.1 GHz, and 
99.3-101.0 GHz. The continuum portion of the Band 6 data contained two spectral windows 
(both within a single baseband), again with dual linear polarizations and eight 250 MHz 
spectral channels across a ~1.8 GHz total band. The frequency range covered by the Band 6 
observations was ~228.7-230.4 GHz. 

Using the continuum emission from the Mira AB binary itself, we performed self- 
calibration on the pipeline data for both ALMA bands. Initially, several iterations of phase- 


^https: / / almascience.nrao.edu/alma-data/science-verification. 
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only self-calibration were performed using a solution interval equal to the record length of 
the data (~6 seconds in Band 3 and seconds in Band 6) until convergence was reached. 
These solutions were applied, and additional iterations were performed, solving for both 
amplitude and phase with a solution interval of 60 seconds. To avoid a possible drift in the 
amplitudes, the gains were normalized during this latter step, and the resulting change in 
the amplitude scales of both bands were measured to be <0.7%. 


3.3. Imaging the Data 

Imaging and deconvolution of the JVLA and the ALMA data were performed using 
IMAGE in AIPS. For the images presented here, we used a Briggs robust value of 77=0, a 
cell size of 5 mas, and a circular restoring beam with a FWHM equal to the mean of the 
dimensions of the dirty beam at the appropriate frequency (see Table 1). Corrections were 
applied during imaging for the frequency dependence of the primary beam and the expected 
spectral index of Mira A within the JVLA band (a=1.86; RM97). Images of the Mira AB 
binary at each of the three observed frequencies are presented in Figures [T], [2l andO 


4. Results 


4.1. The Millimeter Wavelength Flux Densities and Sizes of Mira A and B 

4 . 1 . 1 . Mira A 

We have used two approaches to measure the flux density and size of Mira A based on 
the current data. The hrst was to ht a uniform elliptical disk model to the visibility data 


Table 1. Deconvolved Image Characteristics 


Band 

(1) 

(mas) 

(2) 

0B(min) 

(mas) 

(3) 

PA 

(deg) 

(4) 

(mas) 

(5) 

rms noise 
(/zjy bearn”^) 

(6) 

Q (46 GHz) 

45 

34 

-23 

39.7 

29.0 

3 (96 GHz) 

79 

58 

+42 

68.7 

19.7 

6 (229 GHz) 

32 

23 

+ 19 

27.3 

34.8 


Note. — 0B(maj) and major and minor axis diameters of the 

synthesized dirty beam measured at FWHM using “robust 0” weighting of the 
visibilities. The position angle (PA) of the beam is measured east from north. 
Deconvolved images used for the present analysis used circular restoring beams 
with FWHM as specified in column 5, and where 0^ is equal to the mean of 
0B(maj) and ^^(min). 






Table 2. Measured Photospheric Parameters for Mira A 


Band 

Stellar Phase 

^maj 

^min 

PA 

e 

s„ 

n 



(mas) 

(mas) 

(degrees) 


(mjy) 

(K) 

( 1 ) 

(2) 

(3) 

(4) 

(5) 

(6) 

(7) 

(8) 




Values Based on 

Uniform Elliptical Disk Fits to Visibility Data 



Q (46 GHz) 

0.76 

60.0±2.0 (0.7) 

49.3±1.6 (0.5) 

156±4 (1) 

0.82±0.04 

7.4±1.5 (0.04) 

2110±440 

3 (96 GHz) 

0.47 

53.4±11.2 (0.1) 

42.3±13.0 (0.1) 

142±36 (0.3) 

0.79±0.29 

34.4±6.9 (0.01) 

2900±1200 

6 (229 GHz) 

0.51 

46.2±2.3 (0.03) 

41.5±2.3 (0.02) 

144±11 (0.1) 

0.90±0.07 

150±30 (0.04) 

2620±560 



Values Based 

[ on Elliptical Gaussian Fits to Image Data 



Q (46 GHz) 

0.76 

37.5±2.1 (0.7) 

31.7±2.0 (0.7) 

147±6 (4) 

0.84±0.07 

7.6±1.5 (0.07) 


3 (96 GHz) 

0.47 

32.8±18.3 (0.2) 

26±26(0.2) 

139±3 (1) 

O'Sfoi 

34.4±7.0 (0.04) 


6 (229 GHz) 

0.51 

30.4±0.4 (0.03) 

27.0±0.4 (0.03) 

143.4±0.5 (0.4) 

0.89±0.02 

153±31 (0.09) 



Note. — Stellar phase was derived using the ephemeris available from the database of the American Association of Variable Star Observers 
(AAVSO) at http://www.aavso.org The error bars on the measured quantities include formal, systematic, and calibration uncertainties, 
as described in the Appendix. For measured quantities, the value in parentheses indicates the contribution to the total error budget from 
the formal fit uncertainties. For a uniform disk, ^xnaj s^nd ^min ^^e major and minor axis sizes of the disk, respectively; for a Gaussian 
fit they represent the the FWHM dimensions of the elliptical Gaussian, after deconvolution with the size of the dirty beam. For a resolved 
source, ^maj measured from a Gaussian fit is expected to be 0.625 times the value of ^maj measured from a uniform disk fit. The ellipticity, 
e, is calculated as the ratio of ^min to ^maj ■ The brightness temperature was computed as T 5 = 25i//(fci/'^Tr^niaj ^min) where c is the 
speed of light and k is the Boltzmann constant. 


Table 3. Measured Parameters for Mira B 


Band 

^maj 

^min 

PA 

Su 


(mas) 

(mas) 

(degrees) 

(mJy) 

( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 

Q (46 GHz) 

18±5 (4) 

11±6 (5) 

172±30 (29) 

0.97±0.20 (0.05) 

3 (96 GHz)“ 

20±20 (3) 

19±19 (4) 

46±32 (32) 

2.5±0.5 (0.04) 

6 (229 GHz) 

25.4±0.6 (0.4) 

23.4±0.6 (0.4) 

69±7 (7) 

11.54=2.3 (0.09) 


®'After accounting for systematic errors, Mira B is formally indistinguishable from a 
point source in the Band 3 data. 


Note. — Parameters are derived from elliptical Gaussian fits to the image data de¬ 
scribed in Table 1. The error bars on the measured quantities include formal, systematic, 
and calibration uncertainties, as described in the Appendix. For measured quantities, the 
value in parentheses indicates the contribution to the total error budget from formal fit 
uncertainties. 
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at each observed frequency using the AIPS task OMFIT. The results are summarized in the 
top section of Table 2. The second method was to fit elliptical Gaussians in the image plane 
using the “robust 0” images summarized in Table 1. General expressions for an elliptical 
Gaussian and for the visibility amplitude of a uniform disk may be found in Gondon (1997) 
and Pearson (1999), respectively. For a resolved source, the FWHM size measured from an 
elliptical Gaussian fit (after deconvolution with the size of the dirty beam) is expected to 
be a factor of 0.625 times the major axis of the elliptical disk ht (Pearson 1999), consistent 
with our measurements. The derivation of the error bars given in Table 2 is described in the 
Appendix. 

The sizes derived from our uniform elliptical disk fits are overlaid on the contour images 
presented in Figures [H [21 and El We have resolved the stellar disk at all three of the observed 
frequencies. Furthermore, our measurements point to a trend of decreasing size for the radio 
photosphere with increasing frequency. The arithmetic means of the angular diameters in 
Table 2 translate to radio photospheric radii of 4.5x10^^ cm (3.0 AU) at 46 GHz, 3.9x10^^ cm 
(2.6 AU) at 96 GHz, and 3.6x10^^ cm (2.4 AU) at 229 GHz. The radio photosphere model of 
RM97 predicts that the size of the radio photosphere should decrease at higher frequencies, 
and this is the hrst time that this has been confirmed observationally for an AGB star. A 
similar effect has been previously seen in the red supergiant Betelgeuse (him et ah 1998). 

The angular sizes that we measure for the radio photosphere are intermediate between 
values previously measured for the “molecular layer” of Mira in the infrared (Perrin et al. 
2004) and the sizes of the SiO u=l and u=2, J=l-0 maser emitting “rings”, whose radii 
have been found to range between ~30-40 mas (3.3-4.4 AU) over the course of the stellar 
pulsation cycle (Gotton et al. 2006). These results suggest that the SiO masers arise just 
outside the radio photosphere, consistent with previous hndings (Reid & Menten 2007; see 
also Gray et al. 2009). In addition to resolving the photosphere, our new measurements 
point to a statistically signihcant elongation in the shape of Mira’s radio photosphere at 
a level of ~10-20% (i.e., e ~0.8-0.9), with the major axis of the elongation lying along a 
position angle (PA) between 139 to 156 degrees. 

FigurejUplots the short wavelength radio spectrum of Mira A (and Mira B, discussed be¬ 
low) based on our current measurements. Additional data points at 338 GHz from Ramstedt 
et al. (2014) are also included. The latter measurements were derived from a combination of 
ALMA data taken on 25 February 2014 and 2014 May 3 with a relatively compact antenna 
conhguration that only marginally resolved the two stars. 

The solid line overplotted on Figure S] is not a ht to the data, but rather represents the 
hux density as a function of frequency predicted by the radio photosphere model of RM97: 
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/ p \ 1.86 f D \ ^ 

= (looF^) 

where D is the distance to the star in parsecs. We see that this model provides an excellent 
match to the data for onr adopted distance of 110 pc. 

Based on the parameters derived from our uniform disk hts, we present in column 8 of 
Table 2 the brightness temperature, T^, of Mira computed for each of the three observed fre¬ 
quencies. Although the uncertainties in the derived Tb values are rather large, these values are 
systematically higher than expected from the radio photosphere model of RM97, which pre¬ 
dicts brightness temperatures of RilOOO K. Our values are also higher than the measurement 
of Tb = 1680±250 K previously reported for Mira by Reid & Menten (2007). Temperatures 
higher than those predicted by the RM97 model would imply an increase in opacity owing to 
more efficient collisional ionization of K, Na, and other low-ionization elements. Consistent 
with this, the photospheric radii that we measure are systematically smaller than the RM97 
model predicts (ignoring uncertainties in the adopted stellar distance). 

The brightness temperature measurements in Table 2 represent lower limits to the stellar 
effective temperature of Mira, and indeed, we find that they are comparable to, or smaller 
than most published effective temperature values [e.g., Tgfj ~ 2630 — 2800 K (Haniff et ah 
1995); Teff ~ 2918 — 3192 K (Woodruff et ah 2004)]. However, the Tb values systematically 
exceed the temperature of ~2000 K derived for the molecular layer of Mira by Perrin et ah 
(2004). This is somewhat surprising in light of the expectation that the radio photosphere 
should lie just exterior to this molecular layer (RM97). 


4-1.2. Mira B 

Measurements of the millimeter wavelength flux density for Mira B at each of the three 
observed frequencies are presented in Table 3 and are plotted on the spectrum in Figure 01 
We see from Figure 0] that at 46 GHz and 96 GHz, the measured flux densities of Mira B 
match closely with the relation derived by Matthews & Karovska (2006). However, both of 
the two higher-frequency measurements he significant above the relation. The origin of this 
deviation is unclear. The model developed by Seaquist et ah (1984) and Taylor & Seaquist 
(1984) for radio emission from symbiotic binaries predicts a flattening or downturn in the 
flux density at higher frequencies, opposite to the observed trend. Furthermore, the period 
of the Mira AB binary (~500 yr; Prieur et ah 2002) is too long for orbital modulation to 
have a measurable effect on the spectral index over the course of the one-year period when 
these data were obtained. 
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Based on the Gaussian fits presented in Table 3, it is clear that we have resolved the radio 
surface of Mira B at the highest observed frequency. At 46 GHz, the source also appears 
to be marginally resolved, although the inferred size is only a fraction of the synthesized 
beam and the uncertainties are large. At 96 GHz, the inferred size is consistent with being 
intermediate between the other two bands; however, after consideration of systematic errors, 
Mira B cannot formally be distinguished from a point source at this frequency. 

Based on the ALMA 229 GHz data, we derive a radius for Mira B of Ri 2.0 x 10^^ cm. This 
is many orders of magnitude larger than the expected photospheric radius of a white dwarf. 
This is therefore consistent with the radio emission arising from thermal bremsstrahlung from 
an ionization-bounded hypercompact Hll region surrounding Mira B, embedded in Mira A’s 
neutral wind. However, for a canonical electron temperature of 10,000 K, the measured 
radius is approximately four times larger than the radius that would be required for an 
optically thick blackbody to yield the measured flux density. This, coupled with the apparent 
change in spectral index at high frequencies imply that a more complex model is needed to 
fully describe the observed short-wavelength radio emission. This is not surprising given 
the complex environment surrounding the two stars. In particular, Mira A’s circumstellar 
envelope is highly structured on the scales that are relevant here (e.g., Mayer et al. 2011; 
Ramstedt et al. 2014). 


4.2. Evidence for Non-uniformities on the Surface of Mira A 

The high spatial resolution, high signal-to-noise ratio, and good u-v coverage of the 
three data sets presented here allow the most detailed investigation to date of the brightness 
distribution on the radio surface of Mira A. In Figure Owe present plots of the azimuthally 
averaged visibilities as a function of baseline length for the three observed frequencies. Over¬ 
plotted as thick black lines are the best-htting uniform elliptical disk models as presented in 
Table 2. Models of circular uniform disks are overplotted in red. The solid red line is circular 
model with a disk diameter equal to the major axis of the best-htting elliptical model, while 
the dashed line has a diameter equal to the minor axis of the best-htting elliptical model. 

We see that a simple uniform disk model does not provide a fully adequate characteriza¬ 
tion of the data at any of the three observed frequencies. At 46 GHz and 96 GHz, the uniform 
disk model slightly over-predicts the measured visibility at u-v spacings near 1-2 MA and 
0.5-1 MA, respectively, while under-predicting the values on the longest baselines. Deviations 
from the uniform elliptical disk model are also evident at 229 GHz, with the model htting 
poorly across the longest baselines and the deviations showing a high level of signihcance. 
To see this ehect in the image plane, in Figure [6] we plot the 229 GHz contours for Mira A 
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over an image formed by subtracting the best-fitting elliptical disk model in the visibility 
domain and imaging the residuals. We also produced an analogous image by differencing 
a deconvolved image of the real data with a deconvolved image of the best uniform disk 
model, and the result is virtually indistinguishable, suggesting that the image residuals are 
not solely artifacts of the deconvolution process or limited image hdelity. Contamination 
from a spectral line within our band can also be excluded as a likely cause of the apparent 
brightness non-uniformity, as oxygen-rich red giants typically contain only a handful of weak 
molecular lines over the observed bandwidth (Tenenbaum et ah 2010; Kaminski et ah 2013). 

We attempted to improve the hts to the visibility data using models that include compo¬ 
nents in addition to the uniform disk, including a ring (emulating limb brightening), a point 
source, or a Gaussian. For the 46 GHz and 96 GHz data, we were unable to meaningfully 
constrain models including additional components, and none of these three simple models 
produced a signihcant improvement in the £t (as assessed by visual inspection of the residu¬ 
als and the improvement in the statistic). However, for the 229 GHz data, two-component 
models produced a clear improvement in with the elliptical disk-|-Gaussian ht providing 
the best match to the data (Figure [7]). A ht with a disk plus two point sources was also 
attempted, but was less satisfactory than the elliptical disk-|-Gaussian ht and required an 
additional free parameter. 

For the elliptical disk-|-Gaussian model shown in Figure [TJ the position angle of the 
disk component was hxed at 144.4° (i.e., to the value found originally from htting a disk 
model without additional components; see Table 2). The total hux density in this model 
was 151 mJy, with 39% in a Gaussian component with FWHM 26 mas and a minor-to- 
major axis ratio of 0.92, centered 3.5 mas southwest of the nominal disk center, i.e., at 
the location of the bright region visible in Figure [6l In this model, the major and minor 
axes of the disk component are 13% larger and 2% smaller, respectively, compared with the 
ht given in Table 2. While the elliptical disk-|-Gaussian model is not unique, its success 
nonetheless underscores that the observed visibilities point to a nonuniform, asymmetric 
brightness distribution. 

In its simplest form, the radio photosphere model of RM97 is not expected to exhibit a 
temperature inversion with increasing height, although the dissipation of energy from shocks 
could potentially lead to a region of increased temperature and level of ionization just above 
the radio photosphere. Such a region is expected to become detectable at millimeter and 
submillimeter wavelengths and may manifest itself as detectable limb brightening. While the 
present data do not rule out the presence of a chromosphere-like region, the hts presented 
in Figure [7] suggest that models without azimuthally symmetric limb brightening are able to 
reproduce the data equally well, or better than models that include it. 
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5. Discussion 

We have found excellent agreement between the newly observed millimeter wavelength 
flux densities for Mira A and the radio photosphere model of RM97 (FigureH]), offering strong 
observational support for this model and its extension to higher frequencies. However, we 
note that our current measurement of the flux density of Mira A at 46 GHz is somewhat 
higher than previous measurements at comparable frequencies. Based on 43 GHz data taken 
on 2000 October 26, Reid & Menten (2007) measured S'43 ghz= 4:.8±0.2 mJy, while based 
on data from 2005, Matthews & Karovska (2006) reported values of 2.44±0.44 mJy (on 
May 17) and 2.61±0.37 mJy (on June 28). Matthews & Karovska did not have available 
simultaneous measurements of the SiO maser to use for phase self-calibration (see 03.ip . and 
we estimate that decorrelation losses resulting from the larger residual phase scatter in their 
data may have led to amplitude losses as high as ~ 30%. Nonetheless, even accounting for 
this effect, as well as uncertainties in the absolute flux density scale of up to 20% (typical of 
measurements at this frequency), these data and the earlier measurement of Reid & Menten 
still appear to be consistent with a lower 43 GHz flux density for Mira A compared with our 
latest measurement. 

The measurements of Reid & Menten (2007) and Matthews & Karovska (2006) were 
all obtained near a stellar phase 0 ~0 (i.e., near maximum optical light). Thus while an 
evolution in the radio and millimeter flux density of Mira over timescales of several years 
cannot be excluded, another possibility is that the radio flux density varies over the course 
of the stellar pulsation cycle. Indeed, previous 8.5 GHz measurements by RM97 at several 
phases during the pulsation cycle of Mira implied flux density variations of ~15% at that 
frequency and hinted at a flux density minimum near 0=0. RM97 showed that such a 
trend could be reproduced by a radio photosphere model with low-level shocks (propagation 
speeds ~10 km s“^; see also Reid & Goldston 2002). Additional high-cadence monitoring 
at radio and millimeter wavelengths is clearly needed to settle the question of whether 
Mira’s radio flux density varies appreciably over the course of the stellar pulsation cycle, 
over longer timescales (several years), or both. Further, simultaneous monitoring at several 
frequencies, including longer wavelengths where the stellar disk is unresolved, would help 
to distinguish apparent brightness variations caused by uncertainties in the absolute flux 
density calibration of the JVLA and ALMA at millimeter wavelengths. Galibration at these 
wavelengths, particularly for long baseline measurements, is notoriously difficult since most 
strong, compact radio sources are highly variable at these wavelengths, while planetary 
sources, despite having well determined millimeter flux densities, become resolved on long 
baselines (see RM97 for additional discussion). 

Both the size and ellipticity that we measure for Mira’s radio photosphere at 46 GHz 
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are consistent with the values previously measured at 43 GHz by Reid & Menten (2007), 
although the PA measured by those authors (+39° ± 50°) deviates from our more recent 
value. This hints at the possibility that the shape of the photosphere may change with time. 
However, given the large formal uncertainty in the earlier PA measurement, this difference 
is only marginally signihcant, and additional measurements will be needed to conhrm this 
possibility and to place constraints on the timescale of any such changes. New data sets 
of the quality that we present here would also permit searches for temporal variations as 
small as ~10% in the photospheric radius, even allowing for uncertainties caused by possible 
brightness non-uniformities in the stellar emission. Previous measurements in the near- and 
mid-infrared have already shown evidence of radius variations of this magnitude in Mira over 
the course of the stellar pulsation cycle (~10%; Weiner et ah 2003; Woodruff et ah 2004). 

Measurements at other wavelengths have also supported the existence of an elongation 
of Mira’s photosphere comparable to the one we measure at millimeter wavelengths. For 
example, e ~0.85 has been reported based on optical studies (Karovska et ah 1991; Wilson 
et ah 1992), although those studies suggest that the position angle of the elongation varies 
with time. Using Hubble Space Telescope ultraviolet imaging data, Karovska et ah (1997) 
also reported an elongation along PA?^ 175° (see also Karovska et ah 2005), while Tatebe 
et ah (2008) showed that measurements 11.15/im are consistent with an ellipticity e ~ 0.9 
along PARil69°. 

There are several possible causes for non-circular shapes of radio photospheres, including 
the presence of a close binary companion (see Mauron et ah 2013), non-radial pulsations, 
rotation, or magnetic helds. Our present data do not yet offer strong constraints on these 
various possibilities, but this is expected to change as more extensive databases of radio 
photosphere measurements become available from the JVLA and ALMA in the coming 
years. 

A simple calculation suggests that the tidal force exerted by Mira B would be insufficient 
to cause a distortion of the observed magnitude. One may calculate an order of magnitude 
estimate of the change in radius or “tidal bulge”, AR, raised by the presence of a perturbing 
potential along a line between the center-of-mass of the two bodies as: 


AR 


mi \ s / 


( 2 ) 


where mi is the mass of the primary body, m 2 is the mass of the perturber, R is the radius 
of the star, and s is the separation between the two bodies (e.g., Danby 1992). Assuming 
mi = 2.5 Mq and m 2 = 0.6 Mq (Prieur et al. 2002), R = Ax 10^^ cm fS 14.1.1|) . and s=100 AU 
(Prieur et ah), then Ar 2 x 10® cm. This implies an ellipticity e ~ 5 x 10“®, several 
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orders of magnitude smaller than observed. One possibility is that Mira has an additional 
unseen companion at a much smaller orbital radius. We also note that more sophisticated 
magnetohydrodynamic simulations have shown that a mechanism known as wind Roche- 
Lobe overflow may impose signihcant asymmetries on the shape of Mira’s wind despite the 
large separation between the binary components (Mohamed & Podsiadlowski 2012). In the 
simulations of Mohamed & Podsiadlowski, distortion is seen in the dust shells as close as 
10 AU from the star, although this study does not specihcally address the magnitude of any 
expected distortion on the radio photosphere, which lies interior to the wind launch region. 

Both observations (e.g.. Cotton et al. 2006) and theoretical modeling (Thirumali & 
Heyl 2013) point to magnetic helds as being important in the atmospheric dynamics of Mira 
and in AGB stars in general (see review by Vlemmings 2013). Future contemporaneous 
measurements of radio photospheres and maser emission from various species and transitions 
would be a powerful combination for constraining the roles of magnetic helds and their roles 
in shaping both the photosphere and the stellar outhow over larger scales (see also Gray et 
al. 2009; Amiri et al. 2012). 

The data we have presented here have shown compelling evidence of a nonuniform, asym¬ 
metric brightness distribution on the radio-emitting surface of Mira. Evidence for brightness 
asymmetries on Mira has also been reported previously based on data at a variety of other 
wavelengths (Haniff et al. 1995; Karovska 1999; Stewart et al. 2015). Lim et al. (1998) also 
showed that the red supergiant Betelgeuse {a Ori) exhibits large-scale irregularities in its 
radio surface (see also Richards et al. 2013). A likely cause of these brightness asymmetries 
are large-scale convective cells, which are expected to be present on the surfaces of both 
AGB stars and red supergiants (Schwarzschild 1975; Ghiavassa & Freytag 2015). 


6. Conclusions 

Using data obtained with the JVLA and ALMA, we have presented resolved mea¬ 
surements of the millimeter wavelength continuum emission from both components of the 
Mira AB binary system at 46 GHz (7 mm), 96 GHz (3 mm), and 229 GHz (1 mm). The 
millimeter wavelength flux densities measured for Mira A are consistent with the radio 
photosphere model previously derived by RM97. A comparison with previously published 
measurements of Mira A at 43 GHz suggests that its flux density may vary over the course 
of the stellar pulsation cycle and/or on timescales of years. However, additional measure¬ 
ments are needed to conhrm this possibility. The emission from Mira B is consistent to 
hrst order with free-free emission from a circumstellar hypercompact Hll region of radius 
2 X 10^^ cm, although the millimeter spectrum exhibits an unexpected upturn at higher 
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frequencies, implying that a more complex model is needed to fully describe the emission. 

The stellar disk of Mira A is resolved at all three observed frequencies. We measure a 
mean photospheric radius of ~3.0 AU at 46 GHz and hnd that the radius decreases systemat¬ 
ically with increasing frequency, as expected if the opacity decreases at shorter wavelengths. 
The shape of the radio photosphere of Mira exhibits an elongation, with a flattening of ~10- 
20%. It is presently unclear whether the shape and/or the direction of the elongation evolves 
with time. 

The data presented here reveal evidence of brightness non-uniformities on the radio- 
emitting surface of Mira A. The cause of these non-uniformities may be large-scale convective 
cells, similar to what has been previously proposed to explain the surface irregularities on 
the supergiant Betelgeuse. 

Note added in proof: After this paper was submitted, an independent analysis of the 
ALMA data for Mira AB by Vlemmings et al. (2015) was submitted and accepted for 
publication. The analysis of the ALMA data by these authors differs from our own in that 
for each ALMA band, Vlemmings et al. analyzed the data from the individual observing 
dates separately, and they adopted slightly different flux densities for the ALMA Band 3 flux 
calibrators. They also performed their imaging and data analysis using CASA and a stand¬ 
alone u-v htting program rather than AIPS. The Vlemmings et ah study independently 
conhrms evidence for a hot spot on Mira A and shows that the radio-emitting surface of 
Mira B has now been resolved. However, there are also some notable differences between our 
hndings. First, based on the results of hts to the Band 3 visibility data, there are statistically 
signihcant differences in the size and PA for Mira A between our respective studies; indeed, 
Vlemmings et al. derive a smaller size for Mira A in Band 3 than in Band 6, the opposite 
of what we report here. This in turn leads to a rather different estimate for the Band 3 
brightness temperature. While the respective sizes that we derive for Mira A in Band 6 
agree to within uncertainties, the PA of the major axis differs by more than 90°. We note 
that our value is consistent with the visible elongation of the star in the image plane (e.g.. 
Figure 3). Another important distinction between our results is that Vlemmings et ah report 
a shallower spectral index for Mira A than we have advocated in the present work. They 
also hnd evidence that the spectral index changes signihcantly between 94 and 229 GHz. 
In contrast, we have found a single spectral index of a=1.86 (i.e., the value predicted by 
the radio photosphere model of MR97) to be consistent with recent measurements between 
46 and 338 GHz. Lastly, while both of our studies have revealed evidence for the presence 
of a hot spot on Mira A at approximately the same location, Vlemmings et al. report a 
larger peak residual for this feature after subtracting the best uniform disk £t and infer a 
smaller size for the spot. However, since the inferred size of this feature is only a fraction of 
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the synthesized beam, its size, and hence its inferred brightness temperature, are necessarily 
rather uncertain. 


The authors thank T. Kaminski for valuable discussions and E. Greisen for developments 
in AIPS that aided this work. The JVLA observations presented here were part of NRAO 
program 14A-026. This paper makes use of the following ALMA data: ADS/JAO.ALMA 
#2011.0.00014.SV. ALMA is a partnership of ESO (representing its member states), NSF 
(USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), 
in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by 
ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility 
of the National Science Foundation operated under cooperative agreement by Associated 
Universities, Inc. 


A. Calculation of Uncertainties for Measured Quantities 

The error bars quoted for the measured quantities in Tables 2 and 3 contain contribu¬ 
tions from both formal htting uncertainties and systematic errors, where the total error is 
assumed to be at = (cr| -|- cr^)^'^. The term a/ was taken to be the uncertainties reported by 
the AIPS JMFIT and OMFIT tasks that were used for htting Gaussian and uniform disk 
models, respectively, and which compute uncertainties according to the formalism outlined 
by Condon (1997). The systematic errors, as-, that affect our measurements include terms 
resulting from blurring of the images from phase huctuations on long baselines (i.e., “see¬ 
ing” effects; hereafter ag^p)-, as well as contributions that are difhcult to quantify from hrst 
principles (hereafter (T<j,o), such as those caused by the limited sampling of the u-v plane, by 
the size and shape of the interferometer dirty beam, and by the gridding and deconvolution 
process. 

Assuming that the residual phase noise on long baselines arises from tropospheric phase 
huctuations that follow a Gaussian distribution, then the true size of a source will be ahected 
by convolution with a Gaussian whose FWHM is given by 9p = a'#'8(ln 2), where a = 
a^X/{27rb) (Thompson et al. 1994). Here a^ represents the rms phase huctuation on a 
baseline of projected length b. For the JVLA data (A ~ 7 mm) we hud a^ ~ 29° for the 
three longest baselines (6=34.5 km), implying 9p ~7.5 mas. For the ALMA Band 3 data 
(A ~3 mm), a^ ~ 75° on a baseline with 6=10.66 km, implying 9p ^ 29 mas. Finally, for 
ALMA Band 6 (A Ril.3 mm), a^ ~ 96° for 6=15.2 km, implying 9p f«4.7 mas. Consequently, 
for our size measurements we include in the error budget a relative error term as,p = [9m — 
\/91n ~ ^pM^m where 9m is the measured major or minor axis FWHM. Phase noise also causes 
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a decrease in the correlated amplitude on the longest baselines of magnitude exp(—(T 0 / 2 ) 
(e.g., Thompson et ah). However, the effect on the total integrated flux density of the source 
in our case is quite small, and we have neglected this effect. 

To obtain an estimate of additional contributions as,o to the total systematic error 
budget, we performed a series of numerical experiments. As an initial test, we constructed 
a series of one-dimensional (1-D) Gaussian signals, added to them different levels of random 
noise, and measured the apparent FWHM. We then deconvolved the apparent width from 
the “true” width by subtracting in quadrature a Gaussian representing various observing 
resolutions to yield an intrinsic FWHM. Finally, based on series of 1000 trials for each 
combination of parameters (which bracketed the range of signal-to-noise ratios and source- 
to-beam size ratios encountered in our real data), we computed the rms dispersion in the 
inferred signal width. In this way we have obtained a lower limit to the expected systematic 
uncertainty in the measured source sizes. 

In a second series of experiments, we introduced into our actual visibility data at each 
frequency a series of artihcial elliptical Gaussians and uniform elliptical disks, offset by —0.25 
arcseconds and —1.0 arcseconds south of the phase center, respectively. For the Gaussian 
sources, we set the flux density and axial ratio to the values measured from our hts to Mira A 
(Table 2), but we systematically varied the position angle from —90° to -|-90° degrees, in 
increments of 15°. For the uniform elliptical disks, we performed trials using disk major and 
minor axes equal to our measured values from Table 2, but again varied the position angles 
over several different values. We then measured the parameters of the artihcial Gaussians 
in deconvolved images using JMFIT and of the artihcial disks in the visibility domain using 
OMFIT, respectively, and compared these numbers with the expected values. 

Not surprisingly, we found the errors on the source sizes to be comparable to, or slightly 
higher in our 2-D experiments compared with those derived from our initial 1-D experiments. 
We also found that the 2-D Gaussian hts tended to systematically underestimate the expected 
source hux density and size, typically by a few per cent. This ehect persisted even if smaller 
cell sizes were used to compute the deconvolved images. This ehect may therefore arise 
in part from the intrinsically non-Gaussian shape of the real interferometer beam and/or 
deconvolution biases. 

Based on the above numerical experiments, we adopt as as,o for the image plane and 
visibility domain measurements, respectively, the mean error in each measured quantity from 
all of the 2-D trials at each of the three frequencies. For the sources sizes, the total systematic 
error includes contributions both from as^o and ag^p and is taken to be = (a^p -|- 
In the case of the hux density, we have additionally included contributions to the systematic 
error budget from pointing uncertainties and absolute hux density calibration (assumed to 


19 


contribute relative errors of 1% and 20%, respectively, at all three observed frequencies). 

For Mira A, the uncertainties on the source size and position angle reported in this 
work are dominated by systematic uncertainties, whereas for the weaker radio source Mira B, 
statistical errors dominate at two of the three observed frequencies. For both stars, the un¬ 
certainty in the integrated flux density is dominated by the absolute calibration uncertainty 
inherent at millimeter wavelengths. 
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Fig. 1.— Image of Mira AB obtained with the JVLA at a frequency of 46 GHz on 2014 
February 23. Contour levels are (—1 [absent],l,2,4,8) x290 /iJy beam“^. The lowest contour 
is 10(7. The restoring beam had a FWHM 39.7 mas. The red ellipse indicates the dimensions 
of the best uniform elliptical disk £t to the visibility data for Mira A. 



Fig. 2.— Image of Mira AB obtained from ALMA data at a frequency of 94 GHz. Contour 
levels are (—1 [absent],1,2,4...128) x 197 pJy beam“^. The lowest contour is lOa. The restor¬ 
ing beam had a FWHM of 68.7 mas. The red ellipse indicates the dimensions of the best 
uniform elliptical disk fit to the visibility data for Mira A. 
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Fig. 3.— Image of Mira AB obtained from ALMA data at a frequency of 229 GHz. Con¬ 
tour levels are (—1 [absent],1,2,4---128)x348 /rJy beam“^. The lowest contour is lOa. The 
restoring beam had a FWHM of 27.3 mas. The red ellipse indicates the dimensions of the 
best uniform elliptical disk £t to the visibility data for Mira A. 
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Fig. 4.— Radio spectra for Mira A (solid points) and Mira B (asterisks). The three lower 
freqnency points are from the present stndy and the point at 338 GHz is from Ramstedt et 
al. 2014. Flnx densities were derived from hts to image data (Table 2). The solid line is 
the spectrnm predicted for Mira A by the radio photosphere model of RM97 (Eqnation 2), 
and the dashed line is the spectrnm derived from earlier measnrements of Mira B (between 
8-43 GHz) from the stndy of Matthews & Karovska 2006 (see 
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Fig. 5.— Visibility versus baseline length for the three data sets presented in this paper: 
46 GHz (top); 94 GHz (middle); 229 GHz (bottom). The solid black line shows the best- 
htting uniform elliptical disk model from Table 2. The red lines represent circular uniform 
disk models with diameters equal to the major axes (solid lines) and minor axes (dashed 
lines) of the elliptical disk fits. 
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Fig. 6.— Contour map of the 229 GHz emission from Mira A (reproduced from Figure |3]), 
overplotted on an image of the residuals after subtraction of the best-htting uniform elliptical 
disk model. The residuals range from —1.4 to +0.9 mJy beam“^, with the peak positive 
residual lying slightly southwest of the nominal disk center. In comparison, the peak surface 
brightness of the observed stellar emission is 71.1 mJy beam“^. The restoring beam (FWHM 
27.3 mas) is indicated in the lower left corner. 
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Fig. 7.— Same as Figure IS:, but with additional models overplotted. The black solid line is 
the uniform elliptical disk model from Table 1. The colored lines show models that include a 
uniform elliptical disk plus an additional component: disk+point source (turquoise dash-dot 
line line); disk-|-ring (purple dotted line); disk-|-Gaussian (blue dashed line). 




